Read and display UW-Madison HSRL L1B scanning data#

Preparing your python environment to run this jupyter notebook#

In addition to have JupyterLab installed, you will need to the following Python packages to run this notebook.

  1. xarray

  2. netcdf4

  3. bokeh

  4. holoviews

  5. hvplot

  6. datashader

  7. pyviz_comms

Introduction#

The ground-based University of Wisconsin (UW)-Madison scanning High Spectral Resolution Lidar (HSRL) L1B (similar to a NASA CALIOP L1B product) scanning product consists of the following geophysical 532nm wavelength calibrated measurements:

  1. particulate backscatter [\(\mathrm{m}^{-1}\,\mathrm{sr}^{-1}\)],

  2. particulate linear depolarization ratio,

  3. particulate optical depth,

  4. and attenuated color ratio between the 1064nm and 532nm wavelengths.

The calibrated 1064nm attenuated backscatter [\(\mathrm{m}^{-1}\,\mathrm{sr}^{-1}\)] will be added in the future, along with uncertainties (i.e. standard deviation) of each geophysical variable.

The UW HSRL produces these measurements in three viewing angle configurations.

  • Similar to most ground-based lidar instruments the UW HSRL make continuous “up looking” measurements which we call vertical stare measurements.

  • Next, the UW HSRL is capable of making continuous horizontal stare measurements at an oblique angle between 70° and 90° relative to zenith.

  • Finally, the HSRL makes repeated scanning measurements between 70° and 90°, nominally in 0.5° steps, where each directional scan is conducted within a minute.

These three viewing measurement configurations are organized in their own netcdf4 groups which will be discussed in the subsequent subsections. Here is how you download the sample (prelimary) 1.4 GB L1B netcdf4 file.

import urllib.request
from pathlib import Path

# The remote URL of the sample L1B file
remote_l1b_url_str = (
    "https://ftp.ssec.wisc.edu/pub/incoming/"
    "bagohsrl_20230814T000000_20230815T000000_30.0s_7.5m_1.0deg_L1B_preliminary.nc"
)
Hide code cell source
# Whereto the file will be loaded
l1b_fileP = Path("~/").expanduser() / Path(remote_l1b_url_str).name

# If the file has not been downloaded yet...download the file
if l1b_fileP.exists() is False:
    _, header_obj = urllib.request.urlretrieve(remote_l1b_url_str, str(l1b_fileP))

    # Print how big the file is that you downloaded, just to 
    # make sure that the whole 1.4 GB of data was downloaded
    print(
        "Downloaded {:.1f} GB of data to file {:s}.".format(
            round(int(header_obj["Content-Length"]) / 1024**3, 1),
            str(l1b_fileP)
        )
    )

print(
    f"The downloaded file is available at:\n{l1b_fileP} . "
)
The downloaded file is available at:
/Users/willemm/bagohsrl_20230814T000000_20230815T000000_30.0s_7.5m_1.0deg_L1B_preliminary.nc . 

The L1B file name convention#

The L1B file name consists of the following parts:

<instrument name>_<start date>_<end date>_<time resolution>s_<altitude resolution>m_<angle resolution>deg_L1B<_tag>.nc

For example the specific UW HSRL instrument is designated as the bagohsrl. The file consists of measurements of the 4th of August where the nominal time and altitude resolutions are 30 seconds and 7.5 meters, and the scanning data has a (downsampled) angular resolution of 1 degree. Finally the tag preliminary provides a label of the file.

Geophysical variable naming convention#

Similar to the CALIOP L1B product the geophysical variable names encode the wavelength. The self explanatory variable names are

  • particulate_backscatter_532nm

  • particulate_linear_depolarization_532nm

  • particulate_optical_depth_532nm

  • attenuated_color_ratio_1064nm_532nm

Vertical stare measurements#

Here is how you get the vertical stare data via xarray. The HSRL geophysical variable dimensions are time and mean sea level (MSL) altitude. When you hover your mouse cursor over the variable name under Data variables you will see the full name. The document icon button to the far right, when clicked upon, shows you the attributes of the corresponding variable.

import numpy as np
import xarray as xr

vert_stare_ds = xr.open_dataset(l1b_fileP, group="vertical_stare")
vert_stare_ds
<xarray.Dataset>
Dimensions:                                  (time: 2880, altitude: 2668)
Coordinates:
  * time                                     (time) datetime64[ns] 2023-08-14...
  * altitude                                 (altitude) float64 7.223 ... 1.9...
Data variables:
    particulate_linear_depolarization_532nm  (time, altitude) float64 ...
    optical_depth_532nm                      (time, altitude) float64 ...
    particulate_optical_depth_532nm          (time, altitude) float64 ...
    particulate_backscatter_532nm            (time, altitude) float64 ...
    attenuated_color_ratio_1064nm_532nm      (time, altitude) float64 ...
Attributes:
    documentation:  github.com/ssec/hsrl

Here are the vertical stare geophysical variables plots. The white columns in the images are the time intervals when the HSRL was either making horizontal stares or scanning. The rest of the white areas, above 5km altitude, are low signal to noise ratio (SNR) measurements that were masked out. Note that the particulate optical depth has negative values, since we are in the process of post calibrating the HSRL.

Hide code cell source
import hvplot.xarray
import holoviews as hv

hv.extension("bokeh")
hv.config.image_rtol = 10
hv.opts.defaults(hv.opts.Image(active_tools=["box_zoom"]))
hv.opts.defaults(hv.opts.QuadMesh(active_tools=["box_zoom"]))
hv.opts.defaults(hv.opts.Curve(active_tools=["box_zoom"], show_grid=True))
hv.opts.defaults(hv.opts.Scatter(active_tools=["box_zoom"], show_grid=True))

# Whether to use data shader to aggregate the image pixels when 
# plot the data; refer to 
# https://hvplot.holoviz.org/user_guide/Customization.html#datashading-options
enable_plot_rasterize_bl = True

(
    # Plot the particulate backscatter
    vert_stare_ds.particulate_backscatter_532nm.hvplot(
        x="time", y="altitude", cmap="turbo", clim=(1e-8, 1e-4), logz=True,
        clabel="{:s} [{:s}]".format(
            vert_stare_ds.particulate_backscatter_532nm.attrs["long_name"],
            vert_stare_ds.particulate_backscatter_532nm.attrs["units"]
        ),
        rasterize=enable_plot_rasterize_bl
    ).opts(title=vert_stare_ds.particulate_backscatter_532nm.attrs["long_name"].replace("_", " "))
    
    # Plot the particulate linear depolarization
    + vert_stare_ds.particulate_linear_depolarization_532nm.hvplot(
        x="time", y="altitude", cmap="turbo", clim=(0, 0.6),
        clabel="{:s}".format(
            vert_stare_ds.particulate_linear_depolarization_532nm.attrs["long_name"],
        ),
        rasterize=enable_plot_rasterize_bl
    ).opts(title=vert_stare_ds.particulate_linear_depolarization_532nm.attrs["long_name"].replace("_", " "))
    
    # Plot the particulate optical depth
    + vert_stare_ds.particulate_optical_depth_532nm.hvplot(
        x="time", y="altitude", cmap="turbo", clim=(None, 0.8),
        clabel="{:s}".format(
            vert_stare_ds.particulate_optical_depth_532nm.attrs["long_name"],
        ),
        rasterize=enable_plot_rasterize_bl
    ).opts(title=vert_stare_ds.particulate_optical_depth_532nm.attrs["long_name"].replace("_", " "))
    
    # Plot the attenuated color ratio
    + vert_stare_ds.attenuated_color_ratio_1064nm_532nm.hvplot(
        x="time", y="altitude", cmap="turbo", clim=(0, 1.0),
        clabel="{:s}".format(
            vert_stare_ds.attenuated_color_ratio_1064nm_532nm.attrs["long_name"],
        ),
        rasterize=enable_plot_rasterize_bl
    ).opts(title=vert_stare_ds.attenuated_color_ratio_1064nm_532nm.attrs["long_name"].replace("_", " "))
).opts(title="MAGPIE, Barbados, Ragged Point, 2023-08-14").cols(1)
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

Horizontal stare measurements#

Unlike the vertical stare measurements, the dimensions of the horizontal stare geophysical variable dimensions are time and range, i.e. the optical path distance between the instrument and the atmospheric volume. Within the netcdf4 group two dimensional altitude coordinates are provided, as a function of time.

(Note: The horizontal stare viewing angle will be added in the near future.)

horz_stare_ds = xr.open_dataset(l1b_fileP, group="horizontal_stare")
horz_stare_ds
<xarray.Dataset>
Dimensions:                                  (time: 2880, range: 2668)
Coordinates:
  * time                                     (time) datetime64[ns] 2023-08-14...
  * range                                    (range) float32 7.495 ... 2e+04
    altitude                                 (range, time) float64 ...
    agl_altitude                             (range, time) float64 ...
    altitude_time                            (range, time) datetime64[ns] ...
Data variables:
    particulate_linear_depolarization_532nm  (time, range) float64 ...
    optical_depth_532nm                      (time, range) float64 ...
    particulate_optical_depth_532nm          (time, range) float64 ...
    particulate_backscatter_532nm            (time, range) float64 ...
    attenuated_color_ratio_1064nm_532nm      (time, range) float64 ...
Attributes:
    documentation:  github.com/ssec/hsrl

Here are the horizontal stare geophysical variables plots. The white columns in the images are the time intervals when the HSRL was either making vertical stares or scanning. The rest of the white areas, above 5km altitude, are low signal to noise ratio (SNR) measurements that were masked out. Note that the particulate optical depth has negative values, since we are in the process of post calibrating the HSRL.

Hide code cell source
(
    # Plot the particulate backscatter
    horz_stare_ds.particulate_backscatter_532nm.hvplot(
        x="time", y="range", cmap="turbo", clim=(2e-6, 7e-6), logz=False,
        clabel="{:s} [{:s}]".format(
            horz_stare_ds.particulate_backscatter_532nm.attrs["long_name"],
            horz_stare_ds.particulate_backscatter_532nm.attrs["units"]
        ),
        rasterize=enable_plot_rasterize_bl
    ).opts(title=horz_stare_ds.particulate_backscatter_532nm.attrs["long_name"].replace("_", " "))
    
    # Plot the particulate linear depolarization
    + horz_stare_ds.particulate_linear_depolarization_532nm.hvplot(
        x="time", y="range", cmap="turbo", clim=(0.01, 0.09),
        clabel="{:s}".format(
            horz_stare_ds.particulate_linear_depolarization_532nm.attrs["long_name"],
        ),
        rasterize=enable_plot_rasterize_bl
    ).opts(title=horz_stare_ds.particulate_linear_depolarization_532nm.attrs["long_name"].replace("_", " "))
    
    # Plot the particulate optical depth
    + horz_stare_ds.particulate_optical_depth_532nm.hvplot(
        x="time", y="range", cmap="turbo", clim=(None, 0.8),
        clabel="{:s}".format(
            horz_stare_ds.particulate_optical_depth_532nm.attrs["long_name"],
        ),
        rasterize=enable_plot_rasterize_bl
    ).opts(title=horz_stare_ds.particulate_optical_depth_532nm.attrs["long_name"].replace("_", " "))
    
    # Plot the attenuated color ratio
    + horz_stare_ds.attenuated_color_ratio_1064nm_532nm.hvplot(
        x="time", y="range", cmap="turbo", clim=(0.3, 0.6),
        clabel="{:s}".format(
            horz_stare_ds.attenuated_color_ratio_1064nm_532nm.attrs["long_name"],
        ),
        rasterize=enable_plot_rasterize_bl
    ).opts(title=horz_stare_ds.attenuated_color_ratio_1064nm_532nm.attrs["long_name"].replace("_", " "))
).opts(title="MAGPIE, Barbados, Ragged Point, 2023-08-14").cols(1)

The following plot shows the altitude of the measurements at ranges 5km.

Hide code cell source
altitude_da = horz_stare_ds.altitude
altitude_da.name = "altitude coordinate"

altitude_da.sel(range=5.0e+3, method="nearest").hvplot(x="time", label="range 5km")

Scanning measurements#

The scanning measurement dimensions are angle (i.e. scan angle), range and the scan start time by which each individual scan is indexed. Accompanied with these dimensions are the corresponding altitude and distance coordinates which can be used in plotting routines.

scanning_ds = xr.open_dataset(l1b_fileP, group="scanning")
scanning_ds
<xarray.Dataset>
Dimensions:                                  (start_time: 316, angle: 21,
                                              range: 2668)
Coordinates:
  * start_time                               (start_time) datetime64[ns] 2023...
  * angle                                    (angle) float64 70.0 71.0 ... 90.0
  * range                                    (range) float32 7.495 ... 2e+04
    distance                                 (angle, range) float64 ...
    altitude                                 (angle, range) float64 ...
    agl_altitude                             (angle, range) float64 ...
Data variables:
    particulate_linear_depolarization_532nm  (start_time, angle, range) float64 ...
    optical_depth_532nm                      (start_time, angle, range) float64 ...
    particulate_optical_depth_532nm          (start_time, angle, range) float64 ...
    particulate_backscatter_532nm            (start_time, angle, range) float64 ...
    attenuated_color_ratio_1064nm_532nm      (start_time, angle, range) float64 ...

In the following plots low SNR measurements are masked out.

Hide code cell source
(
    # Plot the particulate backscatter
    scanning_ds.particulate_backscatter_532nm.hvplot.quadmesh(
        x="distance", y="altitude", cmap="turbo", clim=(1.01e-6, 1e-5), logz=True,
        clabel="{:s} [{:s}]".format(
            scanning_ds.particulate_backscatter_532nm.attrs["long_name"],
            scanning_ds.particulate_backscatter_532nm.attrs["units"]
        ),
        rasterize=enable_plot_rasterize_bl,
        xlim=(0, 1.5e+4),
        ylim=(0, 2500)
    ).opts(title=scanning_ds.particulate_backscatter_532nm.attrs["long_name"].replace("_", " "))
    
    # Plot the particulate linear depolarization
    + scanning_ds.particulate_linear_depolarization_532nm.hvplot.quadmesh(
        x="distance", y="altitude", cmap="turbo", clim=(0.0, 0.6),
        clabel="{:s}".format(
            scanning_ds.particulate_linear_depolarization_532nm.attrs["long_name"],
        ),
        rasterize=enable_plot_rasterize_bl,
        xlim=(0, 1.5e+4),
        ylim=(0, 2500)
    ).opts(title=scanning_ds.particulate_linear_depolarization_532nm.attrs["long_name"].replace("_", " "))
    
    # Plot the particulate optical depth
    + scanning_ds.particulate_optical_depth_532nm.hvplot.quadmesh(
        x="distance", y="altitude", cmap="turbo", clim=(None, 0.8),
        clabel="{:s}".format(
            scanning_ds.particulate_optical_depth_532nm.attrs["long_name"],
        ),
        rasterize=enable_plot_rasterize_bl,
        xlim=(0, 1.5e+4),
        ylim=(0, 2500)
    ).opts(title=scanning_ds.particulate_optical_depth_532nm.attrs["long_name"].replace("_", " "))
    
    # Plot the attenuated color ratio
    + scanning_ds.attenuated_color_ratio_1064nm_532nm.hvplot.quadmesh(
        x="distance", y="altitude", cmap="turbo", clim=(0.0, 1.0),
        clabel="{:s}".format(
            scanning_ds.attenuated_color_ratio_1064nm_532nm.attrs["long_name"],
        ),
        rasterize=enable_plot_rasterize_bl,
        xlim=(0, 1.5e+4),
        ylim=(0, 2500)
    ).opts(title=scanning_ds.attenuated_color_ratio_1064nm_532nm.attrs["long_name"].replace("_", " "))
).opts(title="MAGPIE, Barbados, Ragged Point, 2023-08-14").cols(1)

Here is how you make a timeseries plot of the scanning measurements over altitude and distance intervals \((60\mathrm{m}, 100\mathrm{m})\) and \((2980\mathrm{m}, 3020\mathrm{m})\).

Hide code cell source
# The altitude interval
altitude_tpl = (80 - 20, 80 + 20)

# The distance interval
distance_tpl = (3000 - 20, 3000 + 20)

# Compute the mean values of the geophysical variables over the distance and altitude intervals
mean_scanning_ds = \
    scanning_ds.\
    where(
        (distance_tpl[0] < scanning_ds.distance) & (scanning_ds.distance <= distance_tpl[1]), 
        drop=True
    ).\
    where(
        (altitude_tpl[0] < scanning_ds.altitude) & (scanning_ds.altitude <= altitude_tpl[1]), 
        drop=True
    ).\
    mean(dim=["angle", "range"], keep_attrs=True)

# Plot the mean values of the geophysical variables over the distance and altitude intervals
(
    mean_scanning_ds.particulate_backscatter_532nm.hvplot.scatter(
        x="start_time", xlabel="Scan start time", grid=True
    )
    + mean_scanning_ds.particulate_linear_depolarization_532nm.hvplot.scatter(
        x="start_time", xlabel="Scan start time", grid=True
    )
    + mean_scanning_ds.particulate_optical_depth_532nm.hvplot.scatter(
        x="start_time", xlabel="Scan start time", grid=True
    )
    + mean_scanning_ds.attenuated_color_ratio_1064nm_532nm.hvplot.scatter(
        x="start_time", xlabel="Scan start time", grid=True
    )
).opts(title="Scanning measurements over altitude and distance intervals [60, 100]m and [2980, 3020]m").cols(1)